%% Treasury Data Maturities Scatterplot
%more info in readme.txt

clear
close all
clc

%% Add paths
addpath(genpath('Fig_III'))


%% Output folder for figures
Figures_path=strcat('.\Figures');

%% Graphic settings
set(0,'defaultlinelinewidth',2.5);
set(0,'defaultaxesfontsize',14);  %12
set(0,'defaulttextfontsize',12);  %12
set(0,'defaultlinelinewidth',2.5);
rrs = get(0,'screensize');
rrf = get(0,'defaultfigurepos');
defpos = [00, 00, 1200, 500]; %2100 700
set(0,'defaultfigurepos',defpos);

%% Importing data using detectimportoptions
%Get the import options
opts = detectImportOptions('TreasuryData2.csv');
%View which variables we want to load
%disp(opts.VariableNames);
%Choose which variables to load
opts.SelectedVariableNames = {'KYTREASNO', 'MCALDT', 'TMATDT', 'TMTOTOUT'};
opts=setvartype(opts, {'MCALDT', 'TMATDT'}, 'datetime');
opts=setvartype(opts, {'KYTREASNO'}, 'categorical');
opts=setvartype(opts, {'TMTOTOUT'}, 'double');
%load the data
data0 = readtable('TreasuryData2.csv', opts);
%Shift to end-of-month date
data0.MCALDT = eomdate(data0.MCALDT);
data0.timeToMaturity = years(data0.TMATDT - data0.MCALDT);

%% Summarize the data
%create groups by date
[data0.dateGroup, MCALDT] = findgroups(data0.MCALDT);
% computation by groups: splitapply
% workflow chart: https://www.mathworks.com/help/matlab/ref/splitapply.html
meanMat = splitapply(@mean, data0.timeToMaturity, data0.dateGroup);

%replace NaNs
data0.TMTOTOUT(isnan(data0.TMTOTOUT))=0;
tmtotdate = splitapply(@sum, data0.TMTOTOUT, data0.dateGroup);
%Merge and calculate weighted times 
data1=join(data0, table(MCALDT, tmtotdate));
%data1.weights = (data1.TMTOTOUT ./ data1.tmtotdate);
data1.weighted_times = data1.timeToMaturity.*(data1.TMTOTOUT./data1.tmtotdate);
wmeanMat = splitapply(@sum, data1.weighted_times, data1.dateGroup);

%Create output Excel file
OutTable = table(unique(data0.MCALDT), meanMat, wmeanMat);
OutTable.Properties.VariableNames = {'Date', 'MeanTimeToMaturity', 'WtdAvgTimeToMaturity'};
filename = 'patientdata.xlsx';
writetable(OutTable,'Results/MaturityStructureSummary.xlsx')


%% Create binned heatmap
nyears=length(unique(year(data0.MCALDT)));
ntimes = max(ceil(data0.timeToMaturity));

%Need to set custom histogram bins
DATES = datetime(1960,1,1,0,0,0):calyears(1):datetime(2022,1,1,0,0,0);
DATES_plot = datetime(1960,1,31,0,0,0):calyears(1):datetime(2021,01,31,0,0,0);

%Get 2D Histogram
h=histogram2(datenum(data0.MCALDT)+0.001, data0.timeToMaturity, ...
    'XBinEdges',datenum(DATES),...
    'YBinEdges',0:41);
%Histogram counts
COUNTS = h.Values;

%Imagesc from the image processing package is good for plotting counts on a color scale
imagesc(datenum(DATES), 0:40,(COUNTS'))
set(gca,'YDir','normal') 
%Changing settings
colorbar
ax = gca;
ax.ColorScale = 'log';
colormap(flipud(gray));
cb = colorbar; 
cb.Ticks = [0 1 10 50 200 600 1200];
cb.TickLabels = num2cell([0 1 10 50 200 600 1200]);
hold on
plot(datenum(unique(data1.MCALDT)), meanMat, 'g--', "LineWidth", 2)
plot(datenum(unique(data1.MCALDT)), wmeanMat, 'r-', "LineWidth", 2)
datetick('x')
xlim([datenum(datetime(1959, 01, 31)) datenum(datetime(2022, 12, 31))])
ylabel("Years Until Maturity")
legend("Average by Treasury ID", "Value Weighted Average")
%title("Average Maturity of U.S. Treasuries")

ylabel(cb, 'Unique Treasuries')
hYLabel = get(cb,'YLabel');
set(hYLabel,'rotation',-90)
set(gca,'fontsize',16)


%% Print to PDF
L = 26;
W=12;
set(gcf,'PaperUnits','centimeters'); 
set(gcf,'PaperSize',[L W]); 
fig = gcf; 
fig.PaperUnits = 'centimeters';  
fig.PaperPosition = [0 0 L W]; 
fig.Units = 'centimeters'; 
fig.PaperSize=[L W]; 
fig.Units = 'centimeters'; 
savefigure_pdf(strcat(Figures_path,'\Figure_III'));


